%% Filtered Inference Function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [pr_tt1, pr_tt2, pr_tt3, pr_tl1, pr_tl2, pr_tl3, yhatvec, gapvec] = filtered_inference(prmtr,T,y,mstarU,mstarL,vol_bdate,growth_bdate)

% Number of states to keep track of            
           no_st=8;               % max number of periods to be considered; need to adjust lik_fnc to allow more than 8 
           dmnsion=3^no_st;       % max number of cases to be considered 
    
% Parameters, where regimes: 0 normal; 1 L-shaped recession; 2 U-shaped recession 

           p01 = prmtr(1);      %Pr[St=1|St-1=0]
           p02 = prmtr(2);      %Pr[St=2|St-1=0]
           p11 = prmtr(3);      %Pr[St=1|St-1=1]
           p22 = prmtr(4);      %Pr[St=2|St-1=2]      
           mu0 = prmtr(5);      % expansion growth 
           mu1 = prmtr(6);      % L-shaped recession effect 
           mu2 = prmtr(7);      % U-shaped recession effect
           omega1 = 0;           % no bounceback effect for L-shaped recession         
           omega2 = -mu2/mstarU; % bounceback effect for U-shaped recession 	   
           delta = prmtr(8);     % change in trend growth 
           sigma2_pre = prmtr(9)^2; % volatility before break
           sigma2_post = prmtr(10)^2;  % volatility after break
		   
% A MATRIX OF TRANSITION PROBABILITIES 

           PI = [1-p01 - p02, 1-p11, 1-p22;
                    p01, p11, 0;
                    p02, 0, p22];
              
% INITIALIZING THE FILTER WITH STEADY-STATE PROBABILITIES 

           A = [eye(3)-PI; ones(1,3)];
           en=[zeros(3,1);1];
           prob__t = inv(A'*A)*A'*en;  %PR[S_t=0]|PR[S_t=1]|PR[S_t=2],3x1 steady-state probs
                             
                                                                                              
           pr_trf0 = reshape (PI, [numel(PI),1]);
           pr_trf1=[pr_trf0;pr_trf0;pr_trf0];
           pr_trf2=[pr_trf1;pr_trf1;pr_trf1];
           pr_trf3=[pr_trf2;pr_trf2;pr_trf2];
           pr_trf4=[pr_trf3;pr_trf3;pr_trf3];
           pr_trf5=[pr_trf4;pr_trf4;pr_trf4];
           pr_trf=[pr_trf5;pr_trf5;pr_trf5];
                      

           prob__t = reshape ([prob__t, prob__t, prob__t]', [numel([prob__t, prob__t, prob__t]'),1]) .* pr_trf0;
           prob__t = reshape ([prob__t, prob__t, prob__t]', [numel([prob__t, prob__t, prob__t]'),1]) .* pr_trf1;
           prob__t = reshape ([prob__t, prob__t, prob__t]', [numel([prob__t, prob__t, prob__t]'),1]) .* pr_trf2;
           prob__t = reshape ([prob__t, prob__t, prob__t]', [numel([prob__t, prob__t, prob__t]'),1]) .* pr_trf3;
           prob__t = reshape ([prob__t, prob__t, prob__t]', [numel([prob__t, prob__t, prob__t]'),1]) .* pr_trf4;
           prob__t = reshape ([prob__t, prob__t, prob__t]', [numel([prob__t, prob__t, prob__t]'),1]) .* pr_trf5;
           prob__  = reshape ([prob__t, prob__t, prob__t]', [numel([prob__t, prob__t, prob__t]'),1]); % unconditional prob for 3^8 cases
    
            sU = [0;0;1];

            s0 = kron(ones(3^7,1),(kron(sU,ones(3^0,1))));
            s1 = kron(ones(3^6,1),(kron(sU,ones(3^1,1))));
            s2 = kron(ones(3^5,1),(kron(sU,ones(3^2,1))));
            s3 = kron(ones(3^4,1),(kron(sU,ones(3^3,1))));
            s4 = kron(ones(3^3,1),(kron(sU,ones(3^4,1))));
            s5 = kron(ones(3^2,1),(kron(sU,ones(3^5,1))));
            s6 = kron(ones(3^1,1),(kron(sU,ones(3^6,1))));
            s7 = kron(ones(3^0,1),(kron(sU,ones(3^7,1))));

            sUmat = [s7,s6,s5,s4,s3,s2,s1,s0];

            sL = [0;1;0];

            s0 = kron(ones(3^7,1),(kron(sL,ones(3^0,1))));
            s1 = kron(ones(3^6,1),(kron(sL,ones(3^1,1))));
            s2 = kron(ones(3^5,1),(kron(sL,ones(3^2,1))));
            s3 = kron(ones(3^4,1),(kron(sL,ones(3^3,1))));
            s4 = kron(ones(3^3,1),(kron(sL,ones(3^4,1))));
            s5 = kron(ones(3^2,1),(kron(sL,ones(3^5,1))));
            s6 = kron(ones(3^1,1),(kron(sL,ones(3^6,1))));
            s7 = kron(ones(3^0,1),(kron(sL,ones(3^7,1))));

            sLmat = [s7,s6,s5,s4,s3,s2,s1,s0];
 
    Lvec = [omega1, 2*omega1, 3*omega1, 4*omega1, 5*omega1, 6*omega1, 7*omega1];
    Uvec = [omega2, 2*omega2, 3*omega2, 4*omega2, 5*omega2, 6*omega2, 7*omega2];
    Lvec = Lvec(1,1:mstarL);
    Uvec = Uvec(1,1:mstarU);
            
            
    prtt = zeros(T,3);  %storage space for Pr[S_t=0|Y_{t}]
    prtl = zeros(T,3);  %storage space for Pr[S_t=0|Y_{t-1}]
    yhatvec = zeros(T,1);
    gapvec = zeros(T,1);
      
    d = 0; dd=0;
    
        j_iter = 1;
    while j_iter <= T % Loop over time periods
            
            if j_iter > growth_bdate        
                d = 1;    
            end  
            
            if j_iter > vol_bdate        
                dd = 1;    
            end  
            
         
        f_cast1 = y(j_iter)*ones(dmnsion,1) - ones(dmnsion,1)*mu0 - d*ones(dmnsion,1)*delta - sLmat(:,8)*mu1 - omega1*sum(sLmat(:,8-mstarL:7)')' - sUmat(:,8)*mu2 - omega2*sum(sUmat(:,8-mstarU:7)')';
          
        yhat = ones(dmnsion,1)*mu0 + d*ones(dmnsion,1)*delta + sLmat(:,8)*mu1 + omega1*sum(sLmat(:,8-mstarL:7)')' + sUmat(:,8)*mu2 + omega2*sum(sUmat(:,8-mstarU:7)')';
        
        var_L = (1-dd)*sigma2_pre*ones(dmnsion,1) + dd*sigma2_post*ones(dmnsion,1);                                    

        prob_dd=pr_trf .* prob__;
                    % Pr[S_t,S_{t-1},.....,S_{t-7} | I(t-1)]
                    

         tmp=prob_dd;
         tmp=tmp(1:2187)+tmp(2188:4374)+tmp(4375:6561); 
         tmp=tmp(1:729)+tmp(730:1458)+tmp(1459:2187);
         tmp=tmp(1:243)+tmp(244:486)+tmp(487:729);
         tmp=tmp(1:81)+tmp(82:162)+tmp(163:243);
         tmp=tmp(1:27)+tmp(28:54)+tmp(55:81);
         tmp=tmp(1:9)+tmp(10:18)+tmp(19:27);
         tmp=tmp(1:3)+tmp(4:6)+tmp(7:9);

         prtl(j_iter,:)=tmp';          %Pr[S_t|Y_t]
           
             
        pr_vl=(1./sqrt(2.*pi.*var_L)).*exp(-0.5*f_cast1.*f_cast1./var_L).*prob_dd;
                                                                 % 3^8 x 1 Joint density of y_t,S_t,..,S_{t-7} given past information

        pr_val=sum(pr_vl');         % f(y_t|I(t-1)), density of y_t given past information:  This is weighted average of 3^8 conditional densities

        lik=log(pr_val);

        pro_=pr_vl/pr_val; % Pr[S_t,S_{t-1},...,S_{t-7} | I(t-1),y_t]
                      % Updating of prob. with new information, y_t  
                                           
        yhatvec(j_iter,1) = yhat'*pro_;
        gapvec(j_iter,1) = -(sLmat(:,9-mstarL:8)*Lvec' + sUmat(:,9-mstarU:8)*Uvec')'*pro_;              
                      
        
         tmp=pro_;
         tmp=tmp(1:2187)+tmp(2188:4374)+tmp(4375:6561); 
         tmp=tmp(1:729)+tmp(730:1458)+tmp(1459:2187);
         tmp=tmp(1:243)+tmp(244:486)+tmp(487:729);
         tmp=tmp(1:81)+tmp(82:162)+tmp(163:243);
         tmp=tmp(1:27)+tmp(28:54)+tmp(55:81);
         tmp=tmp(1:9)+tmp(10:18)+tmp(19:27);
         tmp=tmp(1:3)+tmp(4:6)+tmp(7:9);

         prtt(j_iter,:)=tmp';           %Pr[S_t|Y_t]


        prob__t=pro_(1:dmnsion/3,1)+pro_(dmnsion/3+1:2*dmnsion/3,1)+pro_(2*dmnsion/3+1:dmnsion,1);
                       % Integrate out S_{t-7}: then you get Pr[S_t, S_{t-1},....,S_{t-6}| Y_t]  
                       %  3^7 x 1

        prob__=reshape([prob__t,prob__t,prob__t]', [numel([prob__t,prob__t,prob__t]'),1]);  

   j_iter = j_iter + 1;
   end % End of loop over time periods         
    
           pr_tt1 = prtt(:,1);
           pr_tt2 = prtt(:,2);
           pr_tt3 = prtt(:,3);

           pr_tl1 = prtl(:,1);
           pr_tl2 = prtl(:,2);
           pr_tl3 = prtl(:,3);
   
   
end % End of filter